Flow, Ordering and Jamming of Sheared Granular Suspensions 
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We study the rheological properties of a granular suspension subject to constant shear stress 
by constant volume molecular dynamics simulations. We derive the system 'flow diagram' in the 
volume fraction/stress plane (<f>,F): at low <j) the flow is disordered, with the viscosity obeying a 
Bagnold-like scaling only at small F and diverging as the jamming point is approached; if the shear 
stress is strong enough, at higher <f> an ordered flow regime is found, the order/disorder transition 
being marked by a sharp drop of the viscosity. A broad jamming region is also observed where, in 
analogy with the glassy region of thermal systems, slow dynamics followed by kinetic arrest occurs 
when the ordering transition is prevented. 
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Under the effects of external drives granular media ex- 
hibit a variety of complex dynamical behaviors. Com- 
paction under shaking is a well studied phenomenon, 
characterized by slow relaxation followed by a jamming 
transition at high volume fraction [l|, 0, Q with deep 
analogies to thermal glasses @, H, U Hi • In addition shear- 
induced transitions from flowing to jammed states 0] or 
to ordered flowing states are observed @, Q ■ In such a 
complex panorama, the rheology of the different phases 
of a system under shear is far from being clarified, and 
even a clear 'flow diagram' locating the regions with or- 
dered flow, disordered flow, and jamming, is missing. 

Within such a perspective, we investigate by molecular 
dynamics (MD) simulations the rheology of a dense gran- 
ular suspension subject to a constant shear stress in a box 
of constant volume. Alikeprevious experiments on non- 
Brownian particles (see [{J [l(| and references therein), 
we opt for a constant shear stress rather that a constant 
shear rate setup (see Q and references therein) because 
jamming may be precluded in the latter case. We derive 
the system 'flow diagram' as obtained by shearing for a 
long time a disordered assembly of grains, ft includes an 
order/disorder transition line, signaled by a sharp drop of 
the effective viscosity rj, as well as a jamming region. Sim- 
ilarly to supercooled liquids, a crossover from a transient 
disordered flow to a stationary ordered one is also found, 
jamming occurring when ordering is kinetically impeded. 
In the regime of disordered flow, the system viscosity, rj, 
is found to grow with the applied shear stress according 
to Bagnold scaling, but it departs from it when the shear 
stress increases. As the packing fraction grows towards 
a critical value, r\ is found to diverge with a power law. 
We also show that, similarly to glassy thermal systems, 
the long time state of the system in the jamming region 
is dependent on the dynamical preparation protocol. 

We consider monodisperse spherical beads of diameter 



d and mass m, initially randomly located between the up- 
per and lower plates of a box of a given size l x X l y X l z . 
The lower plate is immobile and the upper one may move 
along x in response to a given applied shear stress a, 
i.e., under a constant force F — al x l y (periodic bound- 
ary conditions are used along the x and y axes). Grains 
are studied by MD simulations of a well known linear 
spring-dashpot model (L3 of Ref. [ll|), including particle 
rotations and static friction: normal interaction between 
two beads is characterized by the elastic and viscoelas- 
tic constants k n = 2 ■ I0 5 and j n = 50, given in stan- 
dard units [12[ ; static friction is implemented by keeping 
track of the elastic shear displacement throughout the 
lifetime of a contact. The tangential elastic constant is 
k t = 5.7 • I0 4 ; the static friction coefficient /i = 0.1; the 
restitution coefficient is 0.88. The box top and bottom 
plates are made of spherical beads, randomly 'glued' at 
height z — l z — 1/2 + e and at z = 1/2 + e, respec- 
tively, where e is randomly chosen for each bead in the 
set [—1/4, 1/4]. The area fraction of grains on the plates 
is high (large non-physical overlaps are allowed between 
the constituent particles) and their mass is 300m. We 
explored a range of box sizes, and below we illustrate the 
results by considering the case with l x — l y — 16, and 
l z = 8 where finite size effects are absent. The small 
value of l x and l y , taken here to allow simulations on 
time scales long enough, is partially compensated by our 
periodic boundary conditions (l z has a value typical for 
experiments 0, [1, [13]). To avoid gravity-induced com- 
paction, gravitational acceleration of grains is neglected 
as in density-matching liquid experiments [lCl ]. 

For each value of the grain volume fraction <f> and of 
the applied force F we have performed 10 simulations us- 
ing different initial disordered states, always finding the 
same typical velocities as well as the same long time state 
(we have performed simulations up to t — 8000): the sys- 
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FIG. 1: Velocity of the system upper plate, Vghean as a func- 
tion of time, t, for the shown values of the volume fraction, (j>, 
sheared with a constant force F = 4 • 10 4 . The system starts 
from a disordered initial pack. Panel (a) When (j> is small 
enough the system flows in a disordered configuration (see 
text and Fig. [2] left panel) . Panel (b) When <f) is increased 
the disordered flow is only transient, as the system has an 
abrupt transition to an ordered flow with a reduced viscosity 
(see text and Fig. [2] right panel) . Panel (c) At even higher 
4> the transient disordered flow has a transition to a jammed 
configuration where « s hear = 0. 



tern rheology is therefore characterized by the value of F 
and on <j> (or equivalently on grain number N (l9|b In 
the 4> and F range of values we explored, which is close 
to the hard-particle limit as the maximum deformation 
of a particle is Sd/d < 10~ 2 , three qualitatively differ- 
ent regimes are usually observed in the system dynam- 
ics, corresponding to different behaviors of the measured 
upper plate velocity, ^ s hear(</>, F). They are summarized 
in Fig. [T] which shows i> s hear as a function of time, t, at 
increasing values of the volume fraction, = 0.593, 0.608 
and 0.614, for F = 4 • 10 4 . 

When </> is small enough (Fig.QJi), the system flows in a 
stationary disordered state (as the one depicted in the left 
panel of Figj2j) from the initial random configuration and 
^shcar fluctuates around a constant value depending on <f> 
and F. The degree of ordering of the system is usually 
quantified by the amplitude of the peak of the structure 
factor S(k) for, say, k = (0, 0, 27r/rf) ; in the disordered 
region, S(k) has no peaks at all and takes values typical 
to disordered arrangements, as shown in Fig[5]o. 

At higher <p values, when F is strong enough (Fig.[Tp), 
the system is trapped for a certain time in a transient 
state where i> s hear keeps a constant value up to a moment 
when it suddenly jumps to a higher stationary plateau. 
Correspondingly, the system exits the disordered flowing 
state (left panel of Fig[2|) and develops a layered and 
partially ordered structure (shown in the right panel of 
Fig®. 





FIG. 2: Vertical (ZX, top) and horizontal (YX, bottom) 
snapshots of sections of the system for the run of Fig. [T]d. 
The left and the right panels show sections taken at t = 400, 
when the system flow is disordered, and at t = 800, when 
the flow is ordered. The ordered state is characterized by the 
formation of crystal-like layers in the XY plane, shifted along 
z. Dark beads form the boundaries (the shown circles have 
various sizes as they are different cuts of our monodisperse 
spheres) . 



The high <f> scenario is drastically changed if the driv- 
ing force F is not strong enough (see Fig. [T}:): after a 
transient flow the system jams in a state as disordered as 
the initial one, and the upper plate velocity becomes zero. 
The above observations highlight that disordered states 
can be either stationary, or transient. The latter states 
can undergo transitions towards ordered flowing states, 
signalled by a marked increase of the shear velocity, or 
alternatively freeze towards jammed packs. 

As the velocity profiles appear to be approximately 
linear in the considered cases, the system resistance to 
flow can be quantified via an effective viscosity r\ = 
F/v s hcai-(4>: F). In Fig. [5]: we depict the 'flow diagram' of 
the suspension in the (<fi, F)-plane. We first consider the 
flow in the disordered region. Figure [3] shows the depen- 
dence of the effective viscosity on the applied shear force 
F, for several values of the volume fraction of the system, 
4>. At low <fi and small applied shear stress, the viscos- 
ity increases to a good approximation with a power law, 
■q oc F a , where the exponent is a ~ 0.5. Accordingly, 
the relation between shear stress (a oc F) and shear rate 
(7 oc fshoar) is (J oc 7^, with (3 = 1/(1 — a) ~ 2, in 
good agreement with the prediction of Bagnold scaling 
((3 = 2). 

Within the disordered region, for a given force F, the 
dynamics of the system strongly slows down as cj> in- 
creases. This is apparent from Fig. [4] where we show 
that, at constant F, the effective viscosity diverges with 
a power law r){4>) oc (<fi c (F) — <fi)~ b as the packing frac- 
tion increases. Such a relation recorded in the disordered 
regime anticipates the presence of a jamming critical vol- 
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FIG. 3: The effective viscosity, r\ = F/v s hen,r{F), as a func- 
tion of the applied external force F for different values of the 
sample volume fraction (<f) = 0.573 — 0.619), recorded in the 
disordered flow regime. At small F the viscosity increases, at 
a good approximation, as a power law in F a , with a ~ 0.5, 
which implies Bagnold scaling a oc j 13 , with /3 ~ 2. Inset: A 
similar power law scaling is found for rj(F) is the ordered flow 
regime, but almost no dependence with <j> is observed. 
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FIG. 4: At a given force F, the effective viscosity measured 
in the disordered flow regime diverges as a power law as the 
volume fraction of the sample approaches a jamming critical 
threshold 4> C {F). Interestingly, data collected at several F 
values collapse on a master power law, suggesting that its 
exponent is F independent. Here, we have scaled the data 
in such a way that r){(j>c — 0.01) = 1. Inset: the fitted value 
to the critical volume fraction, <f> c , increases with the applied 
external shear force. 



ume fraction 4> C (F) (see Fig. [H inset). Interestingly, the 
fit of 4> C (F) from the small 0- value region always overes- 
timates the volume fraction where we indeed observe the 
system to jam (see Fig. [5]), in analogy to glass forming 
materials [13j. Fig. 0] also shows that data recorded at 
different F collapse one onto the other when plotted as 
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FIG. 5: Panel (a) and (b) show the dependence of the 
effective viscosity 77 and of the structure factor S(k) (for k = 
(0, 0, 2n/d)) on the inverse volume fraction, for F = 

10 5 . In the low <j> regime the system flow is disordered as 
S(k) has typically small values. By increasing (f>, there is 
an initial transient disordered flow (S(k) small) with growing 
77, followed by a sharp transition to an ordered flow (with a 
much higher S(k), see text) with a sharp drop in the viscosity 
(see Fig. [lb) . Panel (c) : The system 'flow diagram' showing 
the long time flow regime as a function of the inverse volume 
fraction and of the applied force F: circles mark the 
region of disordered flow, open squares mark ordered flow, 
and diamonds the jammed region. Stars indicate the inverse 
critical volume fraction cj)" 1 defined from the fit in Fig. [4] 



a function of 4> C {F) — <fi, suggesting that the exponent 
b ~ 0.8 is independent of F. This kind of power law de- 
pendence is frequently found in colloidal systems under 
shear (lq |. with the exponent b varying between 1 and 2. 

We now consider the regime observed at higher cf> 
where, under shearing, the viscosity exhibits a sharp 
transition from a disordered to a faster ordered flow as 
in Fig QJ). The change to the faster flow, with the corre- 
sponding reduction of the effective viscosity, is associated 
to the formation of partially ordered layers parallel to the 
shearing plate (Fig. \2§ , as previously observed in experi- 
ments on granular and colloidal suspensions 8, 14. 151. 
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The connection between ordering and viscosity reduc- 
tion is apparent in Fig. [5ja,,b, where we show the depen- 
dence of both rj and S(k) on the inverse volume fraction 
of the system (squares indicate values recorded in the 
transient regime, circles the asymptotic ones). In the low 
4> region, before the ordering transition, S(k) has small 
values typical to disordered arrangements. At higher </>, 
the system can flow disordered, with increasing viscosity 
and S(k) having still comparable small values. Eventu- 
ally, when a sharp drop in the viscosity is observed, it is 
accompanied by a clear increase of the structure factor 
S(k), pointing out the presence of order in the system, 
as seen in the right panel of Fig. [2] In the ordered flow 
region, rj is almost insensitive to <j) in the range here ex- 
plored (see FigEK) and has an approximate Bagnold-like 
scaling in F, rj cx F a , with a ~ 0.5 (see inset in FigJ3|).In 
the high regime, if the applied shear stress is not strong 
enough, the initial transient disordered flow does not gen- 
erally result in an ordered stationary flow, since a full 
dynamical arrest (see Fig. [5fc) with w s hear = is found. 
Jamming occurs when the system disordered configura- 
tions are kinetically trapped in states where further shear 
(at the given value of F) becomes impossible, the un- 
derlying mechanism being still unclear 17, 3]. Such a 
trapping is broken when F is high enough and ordered 
flow appears. 

The flow diagram of Fig. [SJ; summarizes the properties 
of the long time states reached by the system when 
the initial condition is disordered, and suggests an 
analogy between the (cf>, F) flow diagram and the (<f>, T) 
phase diagram of usual thermal systems, the disordered, 
ordered and jammed states corresponding respectively 
to the liquid, crystalline and glassy phases. To reinforce 
this analogy we have investigated the behavior of a 
system initially prepared in an ordered configuration. 
Interestingly, we have found that in the region of low 
F and high <p, where jamming was previously found, 
the ordered flow is not arrested and keeps going with 
a finite value of rj (although the system might jam at 
times longer than those we can investigate). In the 
other regions of the flow diagram, on the contrary, the 
long time state of the system is that shown in Fig. [SJ; 
regardless of the initial conditions. This dependence on 
the initial conditions is analogous to that observed in 
thermal systems, where the cristalline phase is stable in 
the region where the glassy phase is found. 
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